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We study the synchronous stochastic dynamics of the random field and random bond Ising 
^ chain. For this model the generating functional analysis methods of De Dominicis leads to a 

S , formalism with transfer operators, similar to transfer matrices in equilibrium studies, but with 

dynamical paths of spins and (conjugate) fields as arguments, as opposed to replicated spins. 
In the thermodynamic limit the macroscopic dynamics is captured by the dominant eigenspace 
of the transfer operator, leading to a relative simple and transparent set of equations that are 
easy to solve numerically. Our results are supported excellently by numerical simulations. 
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1. Introduction 



^ \ In spite of the absence of equilibrium phase transitions at finite temperature in 

lO ' one-dimensional Ising chains, the dynamics of such systems (solved formally sev- 

! eral decades ago [ll, 1^) continue to be of interest in the context of ageing phenom- 

ena, see e.g. Disordered versions of such chains, with random bonds and/or 
- random fields, generally require new techniques for solution, unless the disorder 

(<~^ \ can be transformed away as for binary bonds. One method for solving disordered 

chains in equilibrium is based on iteration of partition functions for growing chains, 
constrained by the state of the last spin [3-^. More recentl y su ch models were also 
solved by diagonalisation of replicated transfer matrices [1, 10 1 . The situation with 
■ the dynamics of disordered Ising chains is less satisfactory. Except for special cases, 

j_j \ e.g. [ll|], our analytical methods are still under development, although it is clear 

from numerical simulations (e.g. [l^. [l3l|) and from the equilibrium solution that 
the dynamical phenomenology of disordered Ising chains is rich. A renormalisation 
group approach was advocated in [14]. In this paper we use the generating func- 
tional method of [iH] to handle the disorder, and show that this leads to a transfer 
operator formalism very similar to that found in equilibrium studies; we use paral- 
lel dynamics to keep computations simpler, but results for Glauber dynamics will 



be similar. While developing our study, another study was published [16(], also with 
parallel dynamics but based on the cavity method, which appears to represent an 
alternative but mathematically equivalent perspective on some of our equations. 
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2 A.C.C. Coolen and K. Takeda 

2. Model definitions 

We consider N Ising spins crj G {— 1, 1} on a periodic one-dimensional chain. Their 
dynamics are given by a synchronous stochastic ahgnment to local fields of the 
form hi{a;t) = JjCTj-i + Jj+iO"j+i + ^j(i), with the convention i + N = i for all 
i S {1, . . . , N}, and with cr = (cri, . . . , ajy)- Upon defining pticy) as the probability 
to find the system at time t in state cr, this Markovian process can be written as 

The parameter /3 = > measures the noise in the dynamics, which is fully 
random for /3 = and fully deterministic for /3 — t- oo. The 9i{t) represent external 
fields of the form 9i{t) = 9i + 0i{t), with random frozen parts 6i and weak time 
dependent perturbations 0i{t) that serve to define response functions. The bonds 
Ji and the frozen fields 6i are regarded as quenched disorder, drawn for each site i 
independently from a distribution P{J,6). We write averages over the process ([T]) 
as (■ • ■) and averages over the disorder as Upon removing the time dependent 
parts of the external fields, so that hi{cr; t) = JiOi-i + Jj+iCTj+i + 6*4, the process ([T]) 
obeys detailed balance, and the equilibrium state will be of the Peretto form 

p{a) = z-ie-'^^^C^), (2) 
Hp{a) = -Yl ^^^^ -\Y. log2 cosh[/?/i,(cr)]. (3) 

i i 

The correlation and response functions Cij{t,t') = {ai{t)aj{t')) and Gij{t,t') = 
d{(Ji{t))/d6j{t') will be related by the FDT (fluctuation-dissipation theorem) [l^ 

Gij (r > 0) = -P[Cij (r + 1) - dj (r - 1)] , (r < 0) = 0. (4) 



3. Generating functional analysis 

In order to analyse the macroscopic dynamics of the chain we concentrate on the 
calculation of the disorder averaged generating functional proposed in |15|: 



Z[xl,] = (exp[-i J] 



= E--- E ^k(0),...,^M]expHj] (5) 
cr(o) <T(t„) i t<t„. 

We isolate the local fields at times t € {0, . . . , — 1} in the usual manner via delta 
functions, using the short-hand {dhdh} = Hi [d^i(^)d/ij(t)/27r], which gives 

ZW]= [ {dhdh} Y-.-Yl P('T(0))e^-^[f^>'{'^}] 
cr(o) cr{t„) 

n e''''(*)['''(*)~^''(*)l~*'^'(*)'^'(*)+'^'^'(*+^)'''(*)"^°s2cosh[/3h.(t)]^ j^g-j 
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with the disorder dependent exponent 
J^[{a}, {h}] = ^ logn/ dJdO PiJ, 9)e-'^^ {e/..W+J[/^.{t).._,{t)+/.,_,(iV.(t)]} ^ 

i 

To benefit from the hnear nature of the chain we write ^ in terms of the single- 
site objects <Tj = (fTi(O), . . . ,cTj(tm, - 1)), hj = (/ij(0),... and hj = 
(hj(0), . . . ,hj(tm — 1)), with analogous definitions for di and i/^j. We also introduce 
the time shift matrix S, with entries Sw = dt^t'+i, and the vector u = (l,...,l),so 
that X • Sy = x{t + l)y{t) and u • x = x{t). For factorised and homogeneous 
initial conditions po{cr{0)) = Hi ^'('^i (*-*)) '^^ then obtain 

i <Ti...cr„ 
X JJ(crj,hi,hjlM|<Tj_i,hi_i,hi_i]), (8) 

with a non-symmetric transfer operator M, defined via: 

p(cr(0))e^''-''+^['^-^'^+'^(*™)''(*'"-^W 



(o-,h,h|M|cr',h',h') 



a<ij4^cosh[/3/i(t)]] 
X JdJde P{J,9)e''^''-^-'^^^'^'+^' '^l (9) 



Expression ([8]) is for — )■ c» dominated by the largest eigenvalue Amax of ([9]), 
provided its spectrum is discrete at Amax- In an equilibrium replica analysis [1, 
the relevant kernel would have replicated spins as arguments; here the arguments 
are spin 'paths', field 'paths' and conjugate field 'paths' through time. The fields 
Oi and i/'j were only introduced for generating perturbations, so we may expand 
^[-0] in powers of these fields. To do this efficiently we define 

Tr[i^] = /dhdh J]](cr,h,h|ir|cr,h,h), (10) 
cr 

(cr,h,h|5(t)|cr',h',h') = a{t)6{h-h')5{h-h')6a,c7', (H) 

(cr,h,h|iJ(t)|cr',h',h') = h{t)6{h-h')6{h-h.')6(T,cT' . (12) 



Since Z[0] = 1 for any 6 and Z[0] = Tr[M^] for = 0, according to ^ and ([8]) 
respectively, we can be sure that Tr[M^] = 1 and that all terms of order or order 

~2 

6 in our expansion (which can be written in terms of derivatives with respect to 
6 of .^[0]) must be zero. Thus we may write our expansion in the form 

TTTTT . ■ V- , / ^Tr[5(t)Af^] 1 V- , . X , / f,Tr[S{t)Sit')M^] 
Zm = 1 -iE^«(*)-kpkr -2g^-(^)^-(^ ) Tr[M^] 

-2^1^^^ (*) V'. (t ) ^f^^pr. 

i<j tt' ^ ' 

i^j tt' ^ J 
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We may now use the usual relations to express the quantities of interest in the 
spin chain in terms of derivatives of ([5]), e.g. (crj(t)) = ih™^ dZ[ip]/dijji{t) 
and {ai{t)aj{t')) = - lim^ d"^ Zltj)]/ d^i{t), giving 

m,{t) = {a,{t)) = ^^[^^j \ (14) 
Cij{t,t) = {ai{t)(jj{t')) = Tt[MN] - ^ ' 



a„(t.t') = li,„ ^ = - ,Tr[M^-l.-^5WMl-.lH-(0] 



Left- and right eigenvectors with different eigenvalues of Q are always orthog- 
onal, so we can write ([9]) in the form M = X]a'^^('^)' which the U{X) are 
eigenspace projection operator^, with U{\)U{\ ) = if A 7^ A'. The operator 
= "Y^)^ \^U{X) exchanges dynamical information between sites at distance i. 
Hence |A| < 1 for all A, and since Tr[M^] = 1 for any A^, M must have an eigen- 
value A = 1. Provided the largest eigenvalue is isolated in the spectrum, it follows 
that limAr_s.oo = limAr_^oo X]a ^^U{\) = and the above expressions give 

lim m,(t) = TT[S{t)U{l)]/TT[U{l)l (17) 

TV— >-oo 

lim Cij{t,t') = TT[U{l)S{t)M^-'S{t')]/i:x[U{l)] {i < j), (18) 
lim Gij{t,t') = -iTi[U{l)S{t)M^'-^^H'{t')]/Tr[U{l)]. (19) 

TV— >-oo 

The above quantities are disorder averages of quantities which, by carrying site in- 
dices, will not generally be self-averaging. Hence they will not describe the dynamics 
of an individual realisation of the chain, but averages over many such realisations. 
In contrast, the following quantities are expected to be self-averaging: 

m{t) = IhnN^^j^ Zim{t) =Tr[S{t)U{l)]/Tv[U{l)], (20) 

C{t, t') = limjv-.oo ^ Ei Cii{t, t') = Tv[U{l)S{t)S{t')]/TT[U{l)l (21) 

G(t, t') = lim^v^oo ^ Y.i Gait, t') = -iTr[U{l)S{t)H' {t')]/Tr[U{l)]. (22) 



4. Spectral properties of the transfer operator 

From now on we consider only chains with independently distributed bonds and 
fields, i.e. P{J,9) = P{J)P{9). This is the natural and technically easier scenario. 
To study the spectral properties of M it will be helpful to write this operator as 

(a-,h,h|M|o-:h',h') = {2^)-'-P[<T\hy^-^ jdJdO P(J,0)e-'''"-^-^''[i^'^'+'^' '^l,(23) 
with the probability P[cr|h] of a spin exposed to field path h to follow path cr: 

g/3(7(t+l)/l(t) 

2cosh[/3/i(t)] 



Mt+i)h(t) 

Pk|h]=pKO))n , (24) 



t<t„ 



^If the spectrum of M has continuous parts, the eigenvalue sum becomes an integral. 
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4.1. Reduction of left- and right-eigenvectors 

On the right-eigenvectors of ()23p we carry out the following transformation: 

u^{a,h,h) = I ^^^u;R(a,h,x)e*-P[a|h]. (25) 

Insertion into the eigenvalue equation reveals that WYi{a, h, y) = WB.icr, h — y), and 
after some trivial manipulations we obtain the simplified eigenvalue problem 

XwK{(T,h) = ^ I dh'u;R(cr,'h') jdJdO P{J,e)5[h-eu-JcT']P[(T'\h'+J(T]. (26) 

Writing out the left-eigenvector equation immediately reveals that nL(<T,h, h) = 
UL(<T,h). We now carry out a simple Fourier transformation: 

uU^, h) = l p^^l('t, x)e-'^-. (27) 

Insertion into the left-eigenvalue problem then gives 

XwLic7,h) = ^f dhV(o-',h') jdJde P{J,e)5[h.- Ja']P[a'\h! +e^x+Ja].{2?>) 

The ti;L,R(o',h) represent distributions of field path contributions, conditioned on 
spin paths cr. Given that -P( J, 9) = P{J)P{9) they are connected via 

WKia,h) = Jde Pi9) wU(r,h-0u). (29) 

To see this we simply define the function w{a, h) = JdO P{9) wi,{a, h — 9u) and 
use (f28]l to establish that it obeys 



\w{a,h) = ^f dh!wi^{a',h!) fdJd9d9'PiJ,9)P{9')6[h-e'u-Ja']P[a'\h' +9u+Jct] 

= Y^ f dh' w{a',h') fdJd9 P{J,9)5[h-9u-Ja']P[a'\h' + J(T]. (30) 

Hence w{a, h) obeys (j26p and therefore (j29p holds. We are now left with only one 
eigenvalue problem, and upon combining our results we may summarize 

/■dx e~*'^ 

nL(cr,h,h) = y ^^^^--^(/>(cr,x), (31) 

nR(^,h,h) =P[cr|h] j "^"^^p^^^ j P(^)</'(t,x-0u), (32) 

with (/>((T,x) = WL{cr,x) to be solved from 

A0(cr,h) = /"dh'(^(cr',h') fdJd9 P{J,9)6[h-J(T']P[a'\h'+9u+J(T]. (33) 

For /3 ^ one easily calculates that (/)(<T,h) = 2^*™ Yl,a I '^JP{J)^0^ — J^) and 
that the only possible eigenvalue of ([23|) is A = 1. 
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4.2. Physical meaning of the A = 1 eigenf unctions 

The fields experienced at site i can be writen as hj = + Jj_|_i<Tj_|_i + 0jU, where 
hj^ = JiCTi-i. Apart from the periodicity constraint, all information communicated 
to site i from spins at sites j < i is channeled via h^. The conditional likelihood 
Pj(h^|<T) to observe h?- = at site i, given we know that cTj = cr, thus obeys 

Pi{h^\a) =Y^ fdh' Pi^i{W\a')P[(T'\h' + 0i_iu + Jicr],^(h^- J^cr'). (34) 

CT' 

The spin path cr' at site i — 1 is prescribed in Pj„i(h'|cr'), so Pj_i(h'|cr') no longer 
depends on or J,. Hence if we average over the disorder we obtain 

Pi(hR|cr) = [dW Pi-i{h'\cT') P[a'\h' + 0i^iu + J,a]dih^- Jicr'). (35) 



Disorder averaging removes any site dependence of Pi(h|cr), hence Pj(h|cr) = 
(/)(h|cr), where the latter is now a true conditional probability distribution, al- 
though not corresponding to any specific site, giving the final result 

</.(h|cr)= hm N-^yP,{h\a) = lim N'^Y {6{\i-Jiai_i))\a^=^, (36) 

i i 

(/.(h|cr) = Y [dh'(t>{h'\a') j dJdO P(J, 0)(5(h- Jcr')P[cr'|h' + 0u + Jcr]. (37) 

Equation ([37|) is identical to (p3]) for A = 1. Expression ([36]) obeys causality, i.e. 
(/>(h|(/)) is independent of (T(tmax)- It is reasonable to assume that for A = 1 there 
is only one solution of ()37p that obeys causality, and that non-causal solutions 
will be ruled out by time boundary conditions. Thus we may for A = 1 identify 
</>(cr,h) = (/>(h|cr). Similar arguments underly the cavity approach in (l6| . from 
which p3p can be recovered upon substituting the characteristics of the ID chain. 



5. Calculation of observables 



To calculate the observables (j20|21|22p we need the projection operator U{1). If 
we make the reasonable assumption that for /3 > the A = 1 eigenspace is not 
degenerate, we may use (|31l32p . Tr[C/(l)] = 1, and (/)(cr,h) = 0(h|cr) to write 

(cr, h, h|C/(l)|cr', h', h') = 7~'n^j(cr, h, h)tXL(cr', h', h') 

1 /■ /■ pih(h— X— 6'u)— ih'x' 

= -P[*T|h] j dxdx'0(x|^)0(xV') jdO P{9) ^-p^- , (38) 

with 

7 = ^ / dhdh UR(cr,h,h)uL(cr,h,h) 
cr 

= (27r)-*'" Y [^^ / ^[o-|^u+x+x'](/.(x|cr)0(xV)- (39) 

Since i;^(h|cr) is independent of a{tm) and P[cr|h] is independent of h{tm), we can 
sum in ([39]) over a{tm) and integrate over h{tm) (in that order), resulting in the 
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same expression for 7 but with the replacement tm tm — ^- Further iteration of 
this process leads to 7 = (27r)^*". Hence 

(o-,h,h|f/(l)|o-',h',h') = P[a\h] j dxdx'(^(x|cr)(/.(xV) 

x(27r)-*'" Jd9 P{6y^<^~''-^'''>~'^''''. (40) 

We can now write the dynamical observables (j20l21l22p (using integration by parts 
in the response function, where we take t > t') in the physically transparent form 

m{t) = ^cr(t) / dxdx 4>{x\a)(t>{x'\a) jdO P{e)P[a\eu+x+x'], (41) 

C{t,t') = ^ cr(t)cj(t') / dxdx (/)(x|cr)(/)(xV) [dO P{9)P[a\eu+x+x'], (42) 
a •' •' 

G{t,t') =J2'^i^) y"dxdx' (/>(x|cr)(/)(x» Jde P{e) Jdh 5(h-^u-x-x')^j|^ 
= (3^C{t,t'+l) - Y'^ii) I dxdx (/.(x|cr)(/>(xV) 

X /"d^ P(0)P[crI0u+x+x']tanh[/3(0 + x(O+x'(t'))]}- (43) 

6. Binary bonds and symmetrically distributed random fields 

Our equations take a simpler form when the bonds are binary, i.e. for the choice 
P{J,e) = p(6»)[i(l+ry)5(J-l) + i(l-r?)5(J+l)], where r/G [-1,1]. Insertion into 
([37j) shows that the dynamic order parameter can now be written as 



0(h|o-) = 5]$(^V)J(h-O, (44) 

CT' 

^a'\a) = lim iV"^ V (45) 
i 

with 

$(cr'|cr) = -(! + ??) Yl / Pie)P[a'\a" + eu+(T] 

a" '' 

+ ^(l-r/)^$(-cr"| -cr') [ dO P{e)P[-a'\eu-a" -a]. (46) 
a" 

If, furthermore, we choose random fields with P{—9) = P{6) and unbiased initial 
conditions Po{(t{0)) = \, then /d0P(^)P[-cr'|0u-cr'-cr"] = /d0P(e)P[cr'|eu + 
cr'+cr"], and the operator in (|37|) of which we need eigenfunctions commutes with 
the spin-flip operator (F<^)((t'|<t) = $(-<t'| -cr). We then find that <I>(-cr'[ - cr) = 
$(cr'|cr), and the relatively simple equation 

$((T» = Y ^(o-'V) y d^ P(Q)P{ct'\Bm + o- + (47) 
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C(0,t) 



C{t,t') 




t 

Figure 1. Correlations C{t,t') calculated via numerical solution of II47I I (for t,t' < 10, circles) versus 
correlations measured in numerical simulations with A'^ = 10^ (crosses). Here T = 0.1 and P{J,8) = 
^S(d)[ri5{J—l) + (1— r;)5( J+1)] (which can be mapped onto the non-disordered Ising chain). Left: values of 
C{t, 0), i.e. overlap with the initial state, plotted versus t (which decays as a power law). Right: correlations 
C{t,t') plotted versus the ratio t/t' , for t,t' = 1 . . . 10. Dots show the values of C{t,t') for larger times 
50 < t < t' < 100, as measured in simulations, showing the typical nonequilibrium behaviour C{t,t') ~ 
C{t/t') of the ageing regime. Upper branches: even values of t — t'; lower branches: odd values of t — t' . 



The formulae for the macroscopic observables can now also be simplified. If we use 
the following identity, which follows directly from (I44|47p . 

/ dxdx' 0(x|cr),/.(xV) [dO P(e)P[cr|0u+x+x'] = J] $(o-|o-')«'(cr'|cr), (48) 

then we find m{t) = for all t, and 

C{t, t') = a{t)a{t') <^{a\a'Ma'\a), 
cr <j' 

Git,t')=p[cit,t' + 1)-Y^it) E Hct'\ct)<I> 

a crcT" 

X Ue P{e)P[(T\9u+(T' + a"]tanh[l3{9 + a'{t') + a"{t'))]y (50) 



(49) 



(cT"\a) 



Equations (I47|49|50p no longer depend on rj, since for symmetric field distribu- 
tions and with Jj € {—1,1} the observables C and G are invariant under gauge 
transformations of the type crj — )• Tjcij, so the bonds can be transformed away. 

For P{9) = ^6{d — 9) + ^5{9 + 9), with ^ > 0, we are studying the synchronous 
dynamics random field Ising chain for which the statics was solved in 0]. If 6' > 2 
the dynamics close to T = will be trivial. Each spin freezes into the direction 
dictated by its external field, and the T = order parameter reduces to 



t=i t=i 



(51) 



For 9 < 2 the dynamics remains nontrivial. Since P[<T|h] and ^{cr'\cr) obey causal- 
ity, numerical solution of equations ()47|49|50p is for binary fields quite manageable. 
Examples are shown in figures [T] and [21 and compared to data from simulations 
of chains with = 10^ spins. For those times for which solution of (j47|49|50p is 
feasible, i.e. t,t' < 10, the agreement with simulation data is seen to be perfect. 
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C(0,t) 0, 



C{t,t') 



fS) Jgi 








t 

Figure 2. Correlations C{t,t') calculated via numerical solution of II47I I {t,t' < 10, circles) versus cor- 
relations measured in numerical simulations with N = 10^' (crosses). Here T = 0.1 and P{J,8) = 
l[S(9-^) + 5{e+^)]5(e)[r]S{J-l) + {l-ri)S(J+l)], i.e. weak random fields. Left: values of C(t,0), i.e. overlap 
with the initial state, plotted versus t. Right: correlations plotted versus the ratio t/t' , for t,t' = 1 . . . 10. 
Dots show the values of C{t,t') for larger times 50 < t < t' < 100, as measured in simulations. Upper 
branches: even values of t — t'; lower branches: odd values oit — t'. The system stabilises into a meta-stable 
state on the timescales considered, from which it takes significantly more time to escape. 



7. Beyond single-site observables 



Let us finally turn to the calculation of observables that involve multiple sites. It 
is not difficult to transform the eigenvalue problem for M into an equivalent one 
involving a self-adjoint transfer operator, with the same spectrum. This implies that 
all eigenvalues A of M must be real. We note that, unlike the maximum eigenvalue 
of (p3|) . the non- leading eigenvalues may well depend on the upper time tm-, so 
we will from now on write M{tm) and U{X\tm) instead of M and U{\). Since our 
formulae for observables cannot depend on which choice is made, provided tm is 
equal to or exceeding the largest time argument in the observable, we can combine 
(I17I18I) . the expansion M{t) = XU{X\t), and the property Ti[U{l\t)] = 1. This 
allows us to define and work out the two-site correlation function, with i < j: 



Cij{t,t) = hm \{a,it)a,it)) - {a,{t)) {aj{t))] 



Tr[U{l\t)S{t)M^-'{t)S{t)] -Tv[S{t)U{l\t)] Tt[S {t)U {l\t)] 
X^-\t) TT[U{mS{t)U{X\t)S{t)]. 

A(t)^l 



(52) 



If not only Amax = 1 but also second largest eigenvalue X2(t) is isolated, then the 
evolving correlation length ^(i) in the chain can be calculated via 



1/m = - lim 7 log \n-' J2 Ci,i+L{t, t) 



hm ilog V x'^it) TT[umsmim\t)s{t)] 

L— s>oo L ^ — ' 
A(t)^l 

logA2(t). 



(53) 



So the present transfer operator picture gives a transparent (although not necessar- 
ily trivial) route towards evolving correlation lengths. It is not immediately clear 
whether and how such formulae could be extracted from the cavity formalism (l6| . 
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8. Discussion 



We have shown that apphcation of the elegant generating functional analysis 



method of 15|] to disordered Ising chains (with random fields and/or random 
bonds) leads to a dynamical version of the familiar transfer matrix formalism 
used in equilibrium studies, with a transfer operator whose arguments are spin 
paths, field paths and conjugate field paths. Under weak assumptions (e.g. iso- 
lated largest eigenvalue in the spectrum of the transfer operator) one can take the 
thermodynamic limit and find an exact self-consistency equation for a dynamical 
order parameter, from which disorder-averaged single-site correlation and response 
functions can be calculated explicitly. The latter equation can also be derived from 
cavity arguments [l^, but without the appealing connection with the equilibrium 
transfer matrix formalism. As expected, solving the dynamical order parameter 
equation is still nontrivial. In this paper we have focused on establishing the prin- 
ciples of the method, and we therefore limited ourselves to numerical solution, for 
short times only. In a future study we hope to take further steps, and calculate 
e.g. equilibrium forms for the order parameter (using the parallel dynamics FDT 
relation) as well as the form of the solution in the ageing regime. We have also lim- 
ited ourselves here to investigating those properties than can be extracted from the 
eigenvectors of the operator that correspond to the largest eigenvalue, which implies 
calculating single-site objects only. However, in analogy with the usual procedure 
for equilibrium transfer matrices one can also calculate multi-site quantities (such 
as evolving correlation lengths) from the second largest eigenvalue of the transfer 
operator. Whether and how the same could be done within the cavity formalism of 



16l | is not obvious. Finally, our choice to consider parallel dynamics is not critical 
for the feasibility of the proposed formalism. In the case of sequential (Glauber) 
dynamics one will find a very similar structure, but with transfer operators that 
have continuous time paths rather than discrete time paths as arguments. 
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